공간정보분석
ASS-10
‘1st Law-abiding Dog’
by: 김하정, 서용훈, 황세민
※ 본 문서는 PC 환경에 최적화되어있습니다.
1 서울시 고령자 인구수
1.1 사용한 패키지
1.2 서울시 고령자 인구수 단계구분도
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Maltese-10\data\shp\Seoul_Dong.shp", layer: "Seoul_Dong"
with 424 features
It has 7 fields
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x 179354.7 216408.6
y 436257.1 466546.9
Is projected: TRUE
proj4string :
[+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000
+ellps=GRS80 +units=m +no_defs]
Data attributes:
ADM_DR_CD 동명 Shape_Leng Shape_Area
Length:424 Length:424 Min. : 2323 Min. : 215511
Class :character Class :character 1st Qu.: 3931 1st Qu.: 659043
Mode :character Mode :character Median : 4868 Median : 983122
Mean : 5495 Mean : 1431588
3rd Qu.: 6164 3rd Qu.: 1493250
Max. :20854 Max. :13191610
population objectid 시군구명
Min. : 301.0 Min. : 1.0 Length:424
1st Qu.: 474.2 1st Qu.:106.8 Class :character
Median : 612.0 Median :212.5 Mode :character
Mean : 1536.0 Mean :212.5
3rd Qu.: 796.0 3rd Qu.:318.2
Max. :45455.0 Max. :424.0
old <- read.csv("./data/고령자.csv")
old <- filter(old, !(동 == "소계"))
old <- old[, c(2, 3, 4)]
old2 <- as.character(old$계)
old2 <- gsub(",", "", old2)
old2 <- as.numeric(old2)
old <- cbind(old, old2)
old <- old[, c(1, 2, 4)]
old <- dplyr::rename(old, 동명 = "동")
old3 <- old %>% unite(동명, 자치구, 동명, sep = " ")
map_and_data <- merge(sigungu_shp, old3, by = "동명")
tm_shape(map_and_data) + tm_polygons("old2", id = "동명", convert2density = T, palette = "Blues",
border.col = "black", border.alpha = 0.25)1.3 Moran’s I statistic of spatial autocorrelation
전역적 Moran 1 지수는 특정 공간 단위의 속성 값과 그것의 공간시차와의 상관정도를 나타내 주는 통계치로서 유사한 값들의 전반적인 군집경향을 나타내는 지표이다. 서울시 동별 고령자 수의 공간자기상관을 분석하기 위해 Moran I지수를 살펴보았다.
nb <- poly2nb(map_and_data, queen = TRUE)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)
moran.test(map_and_data$old2, lw, zero.policy = TRUE, na.action = na.omit)
Moran I test under randomisation
data: map_and_data$old2
weights: lw
Moran I statistic standard deviate = 10.698, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.3044727197 -0.0023640662 0.0008226718
공간가중행렬은 Queen contiguity 방식을 사용하였으며, Moran I 지수는 약 0.3044 으로 산출되었다. 또한 모란 산포도를 그려본 결과 양의 관계가 나타나고 있음을 알 수 있다. 따라서 서울시 동별 고령자 인구분포는 공간적 자기상관성을 가지고 있다고 말할 수 있다. 그러나 이 지역안에서도 발생할 수 있는 공간자기상관의 국지적 변이를 고려하여 LISA 중 국지 Moran을 사용하였다.
1.4 local Moran’s I
로컬 모란 지수를 산출한 결과와 그것을 급간으로 나눠 시각화한 결과이다.
local_m <- local_m[, 1]
sigma <- sd(local_m)
mbreak <- c(-Inf, -sigma, 0, sigma, Inf)
group <- cut(local_m, mbreak)
pal <- cm.colors(4)
my.col <- pal[group]
plot(map_and_data, axes = T, col = my.col)1.5 OLS
basic <- read.csv("./data/수급자.csv")
basic <- filter(basic, !(동 == "소계") & !(동 == "기타"))
basic <- dplyr::rename(basic, 동명 = "동")
basic2 <- as.character(basic$가구)
basic2 <- gsub(",", "", basic2)
basic2 <- as.numeric(basic2)
basic <- cbind(basic, basic2)
basic <- basic[, c(2, 3, 15)]
basic10 <- basic %>% unite(동명, 자치구, 동명, sep = " ")
sp <- merge(basic, old, by = "동명")
gf_point(sp$old2 ~ sp$basic2, xlab = "서울시 기초생활수급자", ylab = "서울시 고령자수") %>%
gf_lm(interval = "predict") %>% gf_lm(interval = "confidence")
Call:
lm(formula = sp$old2 ~ sp$basic2)
Coefficients:
(Intercept) sp$basic2
19607.283 8.256
Call:
lm(formula = sp$old2 ~ sp$basic2)
Residuals:
Min 1Q Median 3Q Max
-21582 -5807 -531 4982 33727
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19607.283 662.640 29.590 < 2e-16 ***
sp$basic2 8.256 1.106 7.465 4.73e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8779 on 427 degrees of freedom
Multiple R-squared: 0.1154, Adjusted R-squared: 0.1134
F-statistic: 55.72 on 1 and 427 DF, p-value: 4.727e-13
위에서 시행한 회귀분석의 결과 나온 잔차를 시각화하고 잔차의 공간자기상관을 확인하였다.
bfit <- predict(lm(sp$old2 ~ sp$basic2), data.frame(sp$old2))
resi <- (sp$old2 - bfit)
old3 <- cbind(old3, resi)
old_map <- merge(map_and_data, old3, by = "동명")
tm_shape(old_map) + tm_polygons("resi", id = "동명", palette = "Reds", border.col = "black",
border.alpha = 0.25)
Moran I test under randomisation
data: old_map$resi
weights: lw
Moran I statistic standard deviate = 1.863, p-value = 0.03123
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.0510737063 -0.0023640662 0.0008227506
아주 약한 정도의 자기상관성만을 확인할 수 있다.
1.6 SAR
SAR(spatial autoregressive model)은 종속변수가 공간자기 상관을 갖는 경우로, 공간자기상관의 정보를 독립변수로 회귀식에 설정한다. 현재 회귀모델의 종속변수인 서울시 동별 고령자수의 공간자기 상관은 위에서 확인하였다.
old_map <- merge(old_map, basic10, by = "동명")
mod.sar <- lagsarlm(old_map$old2.x ~ old_map$basic2, data = old_map, listw = lw,
zero.policy = T)
summary(mod.sar)
Call:lagsarlm(formula = old_map$old2.x ~ old_map$basic2, data = old_map,
listw = lw, zero.policy = T)
Residuals:
Min 1Q Median 3Q Max
-20545.7 -5229.1 -1129.1 4417.6 32084.3
Type: lag
Coefficients: (numerical Hessian approximate standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 9334.7590 1339.4500 6.9691 3.190e-12
old_map$basic2 6.6077 1.0180 6.4909 8.534e-11
Rho: 0.45502, LR test value: 61.85, p-value: 3.6637e-15
Approximate (numerical Hessian) standard error: 0.052977
z-value: 8.5891, p-value: < 2.22e-16
Wald statistic: 73.772, p-value: < 2.22e-16
Log likelihood: -4408.147 for lag model
ML residual variance (sigma squared): 63295000, (sigma: 7955.8)
Number of observations: 423
Number of parameters estimated: 4
AIC: 8824.3, (AIC for lm: 8884.1)
res <- mod.sar$residuals
sigma <- sd(res)
classes_fx <- classIntervals(res, n = 5, fixedBreaks = c(-Inf, -sigma, 0, sigma,
Inf), rtimes = 1)
cols <- findColours(classes_fx, pal)
par(mar = rep(0, 4))
plot(map_and_data, col = cols, border = "black")
legend(x = "bottom", legend = names(attr(cols, "table")), cex = 1, fill = attr(cols,
"palette"), bty = "n", ncol = 5, title = "Residuals from SAR Model")1.7 SEM
SEM(spatial error model)은 오류항의 공간자기상관을 통제하기 위해 활용된다. 오류항의 공간자기 상관성은 위에 제시한 바와 같이 매우 낮은 정도에 해당하였으나, SEM을 시행해보았다.
mod.sem <- errorsarlm(old_map$old2.x ~ old_map$basic2, listw = lw, zero.policy = T)
summary(mod.sem)
Call:errorsarlm(formula = old_map$old2.x ~ old_map$basic2, listw = lw,
zero.policy = T)
Residuals:
Min 1Q Median 3Q Max
-21162.4 -5461.5 -1001.6 4550.8 33092.4
Type: error
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) 19397.3185 891.9002 21.7483 < 2.2e-16
old_map$basic2 7.2185 1.1143 6.4779 9.302e-11
Lambda: 0.47357, LR test value: 60.857, p-value: 6.1062e-15
Approximate (numerical Hessian) standard error: 0.05501
z-value: 8.6089, p-value: < 2.22e-16
Wald statistic: 74.113, p-value: < 2.22e-16
Log likelihood: -4408.644 for error model
ML residual variance (sigma squared): 63202000, (sigma: 7950)
Number of observations: 423
Number of parameters estimated: 4
AIC: 8825.3, (AIC for lm: 8884.1)
res <- mod.sem$residuals
sigma <- sd(res)
classes_fx <- classIntervals(res, n = 5, fixedBreaks = c(-Inf, -sigma, 0, sigma,
Inf), rtimes = 1)
cols <- findColours(classes_fx, pal)
par(mar = rep(0, 4))
plot(map_and_data, col = cols, border = "black")
legend(x = "bottom", legend = names(attr(cols, "table")), cex = 1, fill = attr(cols,
"palette"), bty = "n", ncol = 5, title = "Residuals from SEM Model")1.8 AIC
OLS 모형과, SAR, SEM 모형의 적합도를 AIC(Akaike Information Criterion)를 통해 비교하였다.
[1] 9012.198
[1] 8824.294
[1] 8825.287
비교결과, ols 모형보다 SAR 모형과, SEM 모형의 AIC 값이 더 낮은 것을 확인할수 있다.
2 디트로이트의 911 신고전화
주제를 정하던 도중 신고전화 공간분포의 분석을 해보면 좋을 것 같아 해당 주제를 가지고 분석을 해 보았다.
2.1 사용한 데이터셋
사용한 데이터셋으로는 2016년 9월부터 2018년 1월까지 디트로이트의 911 긴급전화 신고 데이터셋과, 디트로이트가 속한 웨인 카운티와 다른 카운티가 같이 섞여있는 센서스 셰이프파일, 그리고 셰이프의 클립연산으로 디트로이트 부분만 추출하기 위한 디트로이트 영역 셰이프 파일이 있다.
모든 자료는 Data Driven Detroit’s Open Data Portal 및 City of Detroit Open Data Portal 에서 구득하였다.
2.2 사용한 패키지
2.3 데이터 전처리
최대한 공간적 자기상관이 드러나게끔 신고전화 수를 지역 인구 수로 나눈 것으로 분석을 진행하였다.
2.3.1 미국 센서스 데이터
디트로이트가 속한 Wayne 카운티 및 Oakland, Macomb 카운티가 혼재된 인구 셰이프를 불러왔다. 또한, 셰이프 객체 하나하나에 얼마나 많은 신고 위치 수가 있는지, 클립연산을 행하기 위해 좌표계 통일이 필요하였다. 따라서 해당 좌표계로 통일하기 위해 해당 좌표계를 추출하였다.
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Maltese-10\data\shp", layer: "KEEDPublishingDemographics"
with 174 features
It has 15 fields
CRS arguments:
+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
2.3.2 911 신고전화 데이터
Municipality(지방자치단체)당 신고건수가 얼마나 있는지 알아보기 위하여 인구 폴리곤과 911 신고위치 데이터를 중첩하여, 폴리곤 객체당 얼마나 많은 수가 있는지 나타내었다.
raw <- read_csv("./data/911_Calls_For_Service.csv")
# raw$call_timestamp <- as.Date(raw$call_timestamp) em2019 <- subset(raw,
# as.Date(raw$call_timestamp) >= as.Date('2018-01-01')) %>% subset(call_timestamp
# <= as.Date('2019-12-31'))
# gcrs <- crs gcrs@projargs <- '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
# #wgs1984
ecpoints <- SpatialPointsDataFrame(raw[, 24:25], raw, proj4string = nad)
# pcsecpts <- spTransform(ecpoints,crs)
spdata@data$counts911 <- poly.counts(ecpoints, spdata)
tmap_mode("view")
tm_shape(spdata) + tm_polygons("counts911", id = "counts911", alpha = 0.6)2.3.3 디트로이트로 지역 국한
폴리곤의 Clip 연산을 수행하여 디트로이트만 나오도록하였다. 구한 디트로이트 시 경계 셰이프는 센서스 셰이프와 CRS를 통일시켜주었다.
bd <- readOGR("./data/shp", layer = "f00e39fd-6902-4d8d-92f7-9f2f5d8156072020329-1-5i4vbz.rjbjg",
encoding = "UTF-8")OGR data source with driver: ESRI Shapefile
Source: "C:\Users\dydgn\OneDrive - 공주대학교\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\Maltese-10\data\shp", layer: "f00e39fd-6902-4d8d-92f7-9f2f5d8156072020329-1-5i4vbz.rjbjg"
with 1 features
It has 11 fields
2.3.4 911 신고전화 - 인구 비
인구당 신고전화 수의 비를 알기 위해 신고 수를 인구 수로 나누었다.
c2p_ratio <- vector(mode = "numeric", length(clip@data$GlobalID))
for (i in 1:length(clip@data$GlobalID)) {
c2p_ratio[i] <- clip@data$counts911[i]/clip@data$Total_Popu[i]
}
c2p_ratio[!is.finite(c2p_ratio)] <- 0
clip@data$c2p_ratio <- c2p_ratio
tm_shape(clip) + tm_polygons("c2p_ratio", id = "NHOOD", alpha = 0.6)2.4 공간적 자기상관
2.4.1 공간가중행렬
공간적 자기상관을 파악하기 위해 필수적인 공간가중행렬을 구축하였다. 구축방식은 Queen’s Case를 사용하였다.
p2n <- poly2nb(clip, queen = TRUE) # 인접행렬
plot(clip, col = "gray", border = "blue", lwd = 2)
plot(p2n, coordinates(clip), col = "red", lwd = 2, add = TRUE)p2n은 폴리곤 데이터 안 인접하는 개체(object)의 목록이다. 예를 들어 첫번째 개체의 인접 목록을 보고 싶다면 다음과 같이…
[1] 9 34
1번 개체의 인접하는 개체 수는 2개가 있다.
[1] 8
[1] 26 129
1번 개체의 부여된 개체코드는 8번이고, 인접하는 두 개체의 코드는 9번과 34번이다.
다음은 인접행렬에 가중치를 부여하여 가중행렬을 만드는 과정이다. 이 사례에서는 인접하는 각각의 개체들에 같은 가중치(1/n)을 부여하였다. 가장 직관적이게 인접 수치를 볼 수 있으나, 연구지역 외곽 개체의 lagged 값에 영향을 주어 잠재적으로 해당 데이터의 진정한 공간자기상관을 overshoot 혹은 undershoot할 수 있다는 단점이 있다.
[[1]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571
2.4.2 Spatially lagged values
인접하는 각각의 개체에 가중치를 할당하였다. 즉, 인접긴급전화비율의 평균을 구할 때, 각각의 개체에 해당 가중치가 적용된다.
2.4.3 Moran’s I (Global)
모란 지수 구하기
복잡한 방법
Lagged된 비율과 원래 비율을 선형회귀에 피팅하여 나온 회귀선의 기울기가 모란 I 지수이므로 위와같이 회귀 모델을 생성하였다.
# plot
plot(x = clip$c2p_ratio, y = call.lag, main = " Moran Scatterplot Call/Pop Ratio")
abline(lm(call.lag ~ clip$c2p_ratio), lty = 3, lwd = 4, col = "red")clip$c2p_ratio
0.2588687
산포도와 모란 I 지수이다.
쉬운 방법
혹은 아래와 같은 함수에 집어넣어서 쉽게 구할 수도 있다.
Moran I test under randomisation
data: clip$c2p_ratio
weights: lw
Moran I statistic standard deviate = 2.8384, p-value = 0.002267
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.258868658 -0.020833333 0.009710448
유의수준 구하기(Monte Carlo Simulation)
복잡한 방법
몬테 카를로 검정에서는 속성값이 임의로 개체에 지정돼서 모란 값을 계산한다. 출력 결과는 속성 값이 연구 영역 전체에 랜덤하게 분포한다는 귀무 가설에서 모란의 I 값을 샘플링한 분포다. 그런 다음 처음에 구한 모란 I 값을 이 표본 분포와 비교한다. 출력 결과는 다음과 같이 나타내어 p-value를 구할 수 있다.
p-value \(= \dfrac{N_{extreme}+1}{N+1}\)
여기서 \(N_{extreme}\)은 처음에 구한 모란 지수보다 분포의 가장자리에 더 가까운 MC출력 모란지수의 수를 의미하며, \(N\)은 MC 시뮬레이션을 몇 번 돌렸는지를 나타낸다.
###### MC
n <- 599 # 실험 수
I.r <- vector(length = n)
for (i in 1:n) {
# 신고전화비율을 임의로 섞기
x <- sample(clip$c2p_ratio, replace = FALSE)
# lag 값을 새로 구하기
x.lag <- lag.listw(lw, x)
# 기울기
M.r <- lm(x.lag ~ x)
I.r[i] <- coef(M.r)[2]
}
# Simulated된 모란 지수 도식 우리가 구한 모란 수 또한 빨간 선으로 도식
hist(I.r, main = NULL, xlab = "Moran's I", las = 1)
abline(v = coef(M)[2], col = "red")[1] 0.01
구한 p 값은 지방자치단체 수준에서 군집화 되지 않을 확률이 적다는 점 (p<0.05)을 시사하며, 귀무가설을 기각한다.
쉬운 방법
혹은 아래와 같이 구할 수도 있다.
Moran I test under randomisation
data: clip$c2p_ratio
weights: lw
Moran I statistic standard deviate = 2.8384, p-value = 0.002267
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.258868658 -0.020833333 0.009710448
Monte-Carlo simulation of Moran I
data: clip$c2p_ratio
weights: lw
number of simulations + 1: 600
statistic = 0.25887, observed rank = 594, p-value = 0.01
alternative hypothesis: greater
몬테카를로 시뮬레이션에서 p-value가 다른 이유는 할 때마다 임의로 섞이기 때문이다. 계속된 니트에서도 p 값은 유의수준인 0.05 밑을 유지하였다.
2.4.4 Moran’s I(Local)
국지적 모란지수의 경우 특정 지역의 값과 인접한 주변 지역들이 갖는 값의 가중 평균 값이 서로 유사하게 나타나면 정적인 자기상관으로, 반대로 특정 지역의 값과 인접한 주변 지역들의 가중 평균값과의 차이가 크게 나타나면 부적인 자기상관으로 나타난다.
국지적 모란지수를 아래와 같이 도식하였다.
clip@data$lom <- localmoran(clip$c2p_ratio, lw)[, 1]
tmap_mode("view")
tm_shape(clip) + tm_polygons(col = "lom", title = "Local Moran's I", midpoint = NA,
legend.format = list(flag = "+")) + tm_style("col_blind") + tm_scale_bar(width = 0.15) +
tm_layout(legend.position = c("left", "bottom"), legend.text.size = 1)2.4.5 LISA 지표
plot(x = scale(clip$c2p_ratio), y = scale(call.lag), main = " Moran Scatterplot 911calls/pop")
abline(h = 0, v = 0)
abline(lm(call.lag ~ clip$c2p_ratio), lty = 3, lwd = 4, col = "red")아래 LISA 지표를 만들어보기 위해 정규화를 한 값을 가지고 산포도를 그려보았다.
2.4.5.1 지표 산정
산포도 상에서 점에 해당하는 값을 1 ~ 4사분면 및 유의수준을 기준으로 반복문 및 다중조건문을 사용하여 분류하였다. 여기서 1은 HH, 2는 LL, 3은 HL, 4는 LH, 5는 유의하지 않음이다.
lom <- localmoran(clip$c2p_ratio, lw)
lomp <- lom[, 5]
s.c2p <- scale(clip$c2p_ratio)
s.lag <- scale(call.lag)
# 반복-조건문
flag <- vector(mode = "numeric", length(clip@data$GlobalID))
for (i in 1:length(clip@data$GlobalID)) {
if (s.c2p[i] >= 0 && s.lag[i] >= 0 && lomp[i] <= 0.05) {
flag[i] <- 1
} else if (s.c2p[i] <= 0 && s.lag[i] <= 0 && lomp[i] <= 0.05) {
flag[i] <- 2
} else if (s.c2p[i] >= 0 && s.lag[i] <= 0 && lomp[i] <= 0.05) {
flag[i] <- 3
} else if (s.c2p[i] <= 0 && s.lag[i] >= 0 && lomp[i] <= 0.05) {
flag[i] <- 4
} else {
flag[i] <- 5
}
}
clip$flag <- flag
2.4.5.2 도식
여기서는 LISA지표를 통해 국지적으로 공간적 연관성을 분석하였다. 인구대비 신고건 수가 많은 다운타운을 중심으로 High-High 클러스터가 나온 것을 볼 수 있다.
breaks <- seq(1, 5, 1)
labels <- c("high-High", "low-Low", "High-Low", "Low-High", "Not Signif.")
np <- findInterval(clip$flag, breaks)
colors <- c("red", "blue", "lightpink", "skyblue2", "white")
plot(clip, col = colors[np])
mtext("Local Moran's I", cex = 1.5, side = 3, line = 1)
legend("topleft", legend = labels, fill = colors, bty = "n")
2.4.5.3 국지적 모란 I p-value
다운타운 주변을 제외한 나머지 지역이 LISA 지표에서 유의하지 않은 이유가 다음의 국지모란지수의 p 값에서 나타난다.
2.5 마치며
신고 건수가 많은 다운타운을 중심으로 클러스터가 형성되어 있으며, 여러 지도 도식을 통하여 핫스팟을 파악할 수 있었다.
하지만, 지역을 디트로이트로 한정하면서 수행한 클립 연산시 중첩한 두개의 폴리곤 사이에서 생긴 Sliver를 처리를 하지 못하여, 가중행렬 생성시 폴리곤 가장자리에 있는 작은 개체들까지 가중치 부여가 되었다.
또한 비율을 구할 때, 몇몇 개체의 인구수 속성 값이 0이라 비율을 구하지 못하고 해당 개체는 비율 값을 0이라 두어서 기울기(모란지수)를 구할 때, 영향을 미쳤다는 한계가 존재한다.
3 서울시 자치구의 Cluster 분석
3.1 사용한 패키지
3.2 자료입력
sp <- readShapePoly("./data/shp/TL_SCCO_SIG_W")
census <- read.csv("./data/Gustats.csv", header = T)
geocensus <- merge(sp, census, by = "SIG_KOR_NM")
xa = coordinates(sp)[, 1]
ya = coordinates(sp)[, 2]
dis = as.matrix(dist(cbind(xa, ya)))
disi = 1/dis
diag(disi) = 0
r1 = census$alien/census$pop # 외국인비율
r2 = census$pop/census$hh # 세대당인구
r3 = census$old/census$pop # 고령인구비율3.3 비계층 군집분석(pam) 결과
서울의 25개 자치구 중에서 외국인비율과 고령인구비율을 토대로 5개 그룹으로 군집분석한 결과는 다음과 같다.
geocensus$r1 = geocensus$alien/geocensus$pop
geocensus$r2 = geocensus$pop/geocensus$hh
geocensus$r3 = geocensus$old/geocensus$pop
pc = geocensus[, c("SIG_KOR_NM", "r1", "r3")]
qq = pam(cbind(geocensus$r1, geocensus$r3), k = 5)
pc$rank1 = order(pc$r1)
pc$rank3 = order(pc$r3)
pc$cluster = qq$clustering
plot(pc$r3, pc$r1, type = "n", xlab = "고령인구비율", ylab = "외국인비율", main = "서울 자치구 인구속성")
text(pc$r3, pc$r1, labels = pc$SIG_KOR_NM, cex = 0.6, pos = 3)
points(pc$r3[pc$cluster == "1"], pc$r1[pc$cluster == "1"], pch = 15, col = "red")
points(pc$r3[pc$cluster == "2"], pc$r1[pc$cluster == "2"], pch = 15, col = "blue")
points(pc$r3[pc$cluster == "3"], pc$r1[pc$cluster == "3"], pch = 15, col = "green")
points(pc$r3[pc$cluster == "4"], pc$r1[pc$cluster == "4"], pch = 15, col = "grey")
points(pc$r3[pc$cluster == "5"], pc$r1[pc$cluster == "5"], pch = 15, col = "violet")
legend("topright", c("Cluster no.1", "Cluster no.2", "Cluster no.3", "Cluster no.4",
"Cluster no.5"), lty = 1, lwd = 2, col = c("red", "blue", "green", "grey", "violet"),
text.col = c("red", "blue", "green", "grey", "violet"))- 제1그룹 : 강북구, 도봉구, 은평구, 중랑구, 노원구
- 제2그룹 : 중구, 종로구, 동대문구
- 제3그룹 : 광진구, 마포구, 성동구, 관악구, 동작구, 성북구
- 제4그룹 : 영등포구, 금천구, 구로구, 용산구
- 제5그룹 : 송파구, 강남구, 서초구, 양천구, 강서구, 강동구
4 서울시 자치구의 공간적 자기상관 분석
4.1 x, y 좌표에 의한 Moran’s I 산정
$observed
[1] 0.03064059
$expected
[1] -0.04166667
$sd
[1] 0.03281405
$p.value
[1] 0.02755632
$observed
[1] -0.009935369
$expected
[1] -0.04166667
$sd
[1] 0.03138569
$p.value
[1] 0.3120109
$observed
[1] 0.1088674
$expected
[1] -0.04166667
$sd
[1] 0.03227261
$p.value
[1] 3.094388e-06
RESULT : 서울 자치구의 x, y 좌표에 의하여 산정한 거리를 토대로, Moran’s I 산정결과는 다음과 같다. 외국인비율과 고령인구비율은 공간적 자기상관이 존재한다. 세대당인구는 공간적 자기상관이 없다.
4.2 삼각망 자료에 의한 Moran’s I 산정
4.2.1 Delaunary 삼각망에 의한 공간적 자기상관
coords = cbind(xa, ya)
ind <- sapply(slot(sp, "polygons"), function(x) slot(x, "ID"))
trinb <- tri2nb(coords)
plot(sp, border = "white", col = "green")
plot(trinb, coords, add = TRUE)
title(main = "Delaunay Triangulation")
text(xa, ya, pos = 3, labels = census$SIG_KOR_NM, cex = 0.6, col = "red")
Moran I test under randomisation
data: r1
weights: nb2listw(trinb, style = "W")
Moran I statistic standard deviate = 2.6625, p-value = 0.003879
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.26640488 -0.04166667 0.01338862
Moran I test under randomisation
data: r2
weights: nb2listw(trinb, style = "W")
Moran I statistic standard deviate = 0.2804, p-value = 0.3896
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-0.01067300 -0.04166667 0.01221797
Moran I test under randomisation
data: r3
weights: nb2listw(trinb, style = "W")
Moran I statistic standard deviate = 4.6466, p-value = 1.687e-06
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.48688034 -0.04166667 0.01293873
RESULT : 서울 자치구의 Delaunary Triangular Link에 의하여 산정한 Neighbor를 토대로, Moran’s I를 산정한 결과는 다음과 같다. 외국인비율과 고령인구비율은 공간적 자기상관이 존재한다. 세대당인구는 공간적 자기상관이 없다.
4.2.2 KNN(N=4) 삼각망에 의한 공간적 자기상관
knn4 = knearneigh(coords, k = 4)
knn4.nb = knn2nb(knn4)
plot(sp, border = "white", col = "green")
plot(knn4.nb, coords, add = TRUE)
title(main = "KNN : K Nearest Neighbours \n (k=4) Triangular Links")
text(xa, ya, pos = 3, labels = census$SIG_KOR_NM, cex = 0.6, col = "red")
Moran I test under randomisation
data: r1
weights: nb2listw(knn4.nb, style = "W")
Moran I statistic standard deviate = 1.9894, p-value = 0.02333
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.19996175 -0.04166667 0.01475247
Moran I test under randomisation
data: r2
weights: nb2listw(knn4.nb, style = "W")
Moran I statistic standard deviate = 0.75438, p-value = 0.2253
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.04591695 -0.04166667 0.01347921
Moran I test under randomisation
data: r3
weights: nb2listw(knn4.nb, style = "W")
Moran I statistic standard deviate = 4.0471, p-value = 2.592e-05
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.44167741 -0.04166667 0.01426315
RESULT : 서울 자치구의 좌표에서 KNN(k=4)를 적용하여 산정한 Neighbor를 토대로, Moran’s I를 산정한 결과는 다음과 같다. 외국인비율과 고령인구비율은 공간적 자기상관이 존재한다. 세대당인구는 공간적 자기상관이 없다.
이후부터는 고령인구비율을 집중적으로 다루뤄보겠다.
4.3 삼각망 자료에 의한 Local Moran’s I 및 G-통계량 산정
4.3.1 Delaunary 삼각망을 이용한 Local Moran’s I 산정
nb1 <- tri2nb(coords)
lom <- localmoran(r3, nb2listw(nb1, style = "W"))
moran.plot(r3, nb2listw(nb1, style = "W"), pch = 19, xlab = "고령인구비율", labels = as.character(sp$SIG_KOR_NM),
main = "Moran Scatter Plot")moran.plot(as.vector(scale(r3)), nb2listw(nb1), pch = 19, xlab = "고령인구비율",
labels = as.character(sp$SIG_KOR_NM), main = "Moran Scatter Plot as Vector") Ii E.Ii Var.Ii Z.Ii
Min. :-0.20772 Min. :-0.04167 Min. :0.1007 Min. :-0.3675
1st Qu.: 0.03169 1st Qu.:-0.04167 1st Qu.:0.1237 1st Qu.: 0.1858
Median : 0.19441 Median :-0.04167 Median :0.1559 Median : 0.6713
Mean : 0.48688 Mean :-0.04167 Mean :0.1653 Mean : 1.2745
3rd Qu.: 0.60626 3rd Qu.:-0.04167 3rd Qu.:0.2042 3rd Qu.: 2.0422
Max. : 2.00084 Max. :-0.04167 Max. :0.2847 Max. : 5.1735
Pr(z > 0)
Min. :0.0000001
1st Qu.:0.0205659
Median :0.2510101
Mean :0.2462355
3rd Qu.:0.4263032
Max. :0.6433708
1상한의 강북구는 주변지역과 함꼐 Local Moran’s I 가 높은 값을 나타내는 지역으로 판단된다.
4.3.1.1 Local Moran’s I 분석
pr = data.frame(cbind(lom[, 1], (r3 - mean(r3))/sd(r3)))
plot(pr, xlab = "Moran's I", ylab = "Standardised Rate", xlim = c(-0.5, 2.5), pch = 19)
text(pr[, 1], pr[, 2], pos = 1, sp$SIG_KOR_NM, cex = 0.6)
abline(h = 0, col = "red")
abline(v = 0, col = "red")
title("Plot : Local Moran's I")
mtext(side = 3, "- 서울 자치구 고령인구비율 -")Local Moran’s I 산정결과는 전체적으로 높은 값을 보이고 있는 것으로 파악된다.
4.3.1.2 Local Moran’s I 그래프
ll = localmoran(r3, nb2listw(nb1, style = "W"))
pa = c("#ffffff", "#f0ff00", "#ffce00", "#ff9a00", "#ff0000")
qk = classIntervals(round(ll, 4), n = 5)
qkp = findColours(qk, pa)
plot(sp, col = qkp)
title("Local Moran's I")
legend("topleft", fill = attr(qkp, "palette"), cex = 0.8, legend = names(attr(qkp,
"table")), bty = "n")
text(coords, labels = sp$SIG_KOR_NM, cex = 0.6)
mtext(side = 3, "- 고령인구 비율 -")4.3.2 KNN 삼각망을 이용한 Local Moran’s I 그래프
비계층적 군집분석 방법인 K-Means(K=5)에 의한 Local Moran’s I에 따른 분포패턴은 다음과 같다.
nb2 <- knn2nb(knn4)
zz = (r3 - mean(r3))/sd(r3)
qq = nb2listw(nb2)
wz = lag.listw(qq, zz)
mlm = lm(wz ~ zz)
aa = mlm$coef[1]
bb = mlm$coef[2]
plot(zz, wz, xlab = "Rate", ylab = "Spatial Lag of Rate")
abline(aa, bb)
abline(h = 0, lty = 2)
abline(v = 0, lty = 2)
title(paste("Moran Plot I =", format(round(bb, 4))))
text(zz, wz, pos = 3, sp$SIG_KOR_NM, cex = 0.6)1상한(주변지역과 같이 높은 값) : 강북구, 도봉구, 중구, 종로구, 은평구, 서대문구, 용산구, 성북구, 노원구
2상한(주변지역에 비해 낮은 값) : 해당 없음
3상한(주변지역과 함께 낮은 값) : 강남구, 송파구, 서초구, 광진구, 양천구, 마포구, 강서구, 가동구, 성동구, 영등포구, 관악구, 금천구, 구로구
4상한(주변지역에 비해 높은 값) : 동대문구, 중랑구, 동작구
4.3.3 Local G 통계량 산정
lgg = localG(r3, nb2listw(nb1, style = "W"))
pa = c("#ffffff", "#f0ff00", "#ffce00", "#ff9a00", "#ff0000")
qk = classIntervals(round(lgg, 4), n = 5, style = "kmeans")
qkp = findColours(qk, pa)
plot(qk, pal = pa, xlab = "고령인구", ylab = "누적확률", main = "누적확률 Plot (K-Means)")plot(sp, col = qkp)
title("Local G-Statistics : K-Means")
legend("topleft", fill = attr(qkp, "palette"), cex = 0.8, legend = names(attr(qkp,
"table")), bty = "n")
text(coords, labels = sp$SIG_KOR_NM, cex = 0.6)
mtext(side = 3, "- 고령인구 비율 -")한강 이남 vs 한강 이북, 구도심 vs 신도심 간의 명확한 차이가 관찰된다
4.3.4 Local Moran’s I 및 Local G 산정결과 출력
Local Moran’s I 및 Local G 통계량을 출력한 결과는 다음과 같다. 또한 Local Moran’s I의 유의확률이 \(\alpha\)=0.05 이하를 추출하면 도봉구, 은평구, 종로구, 강북구, 강남구, 서초구, 송파구 등 7개 지역이 추출되었다.
aaa = data.frame(round(cbind(lom, lgg), 3), row.names = sp$SIG_KOR_NM)
colnames(aaa) <- c("Moran-I", "E(I)", "V(I)", "Z(I)", "Pr(I)", "G-st.")
print(aaa) Moran-I E(I) V(I) Z(I) Pr(I) G-st.
도봉구 1.838 -0.042 0.285 3.523 0.000 2.303
은평구 0.682 -0.042 0.124 2.057 0.020 2.017
동대문구 -0.118 -0.042 0.156 -0.193 0.577 -0.312
동작구 -0.022 -0.042 0.156 0.049 0.481 -0.874
금천구 0.032 -0.042 0.156 0.186 0.426 -1.012
구로구 0.020 -0.042 0.204 0.138 0.445 -1.213
종로구 1.288 -0.042 0.156 3.368 0.000 2.854
강북구 2.001 -0.042 0.156 5.173 0.000 2.539
중랑구 -0.208 -0.042 0.204 -0.367 0.643 -0.517
강남구 1.558 -0.042 0.204 3.540 0.000 -2.588
강서구 0.184 -0.042 0.204 0.500 0.308 -0.681
중구 0.437 -0.042 0.156 1.211 0.113 1.009
[ reached 'max' / getOption("max.print") -- omitted 13 rows ]
Moran-I E(I) V(I) Z(I) Pr(I) G-st.
도봉구 1.838 -0.042 0.285 3.523 0.000 2.303
은평구 0.682 -0.042 0.124 2.057 0.020 2.017
종로구 1.288 -0.042 0.156 3.368 0.000 2.854
강북구 2.001 -0.042 0.156 5.173 0.000 2.539
강남구 1.558 -0.042 0.204 3.540 0.000 -2.588
서초구 0.606 -0.042 0.101 2.042 0.021 -1.670
송파구 1.590 -0.042 0.204 3.611 0.000 -2.590